{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "Created on Tue Dec 01 21:05:42 2015\n",
    "\n",
    "@author: Tim\n",
    "\"\"\"\n",
    "import matplotlib.pyplot as P\n",
    "from pylab import *                ## import scientific database\n",
    "close(\"all\")                       ## close all windows\n",
    "from netCDF4 import Dataset\n",
    "import numpy as np\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "from matplotlib.font_manager import FontProperties\n",
    "from wrf import getvar, interplevel, to_np, get_basemap, latlon_coords, CoordPair, vertcross\n",
    "import glob\n",
    "import cmocean\n",
    "import cmocean.cm as cmo\n",
    "from functions_import import import_VNP02, import_cldprop_l2, import_cldprop_uc_l2, import_MOD021KM, import_MOD06_L2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BRING IN SATELLITE DATA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Case\n",
    "case2 = \"0313\"\n",
    "\n",
    "if case2 == \"0328\" or case2 == \"0329\":\n",
    "    case = \"0328\"\n",
    "else:\n",
    "    case = \"0313\"\n",
    "\n",
    "### Head directory\n",
    "head = \"/glade/derecho/scratch/tjuliano/doe_comble/satellite/data/\" + case + \"/\"\n",
    "\n",
    "### Satellite file names, --> Imagery, Imagery geo location, cloud retrival, cloud retrival geo location, folder name ###\n",
    "#sats = [[\"VJ102IMG\",\"VJ103IMG\",\"CLDPROP_L2_VIIRS_NOAA20\", \"VJ103MOD\", \"viirs_temp\"], # NOAA-20\n",
    "#        [\"VNP02IMG\",\"VNP03IMG\",\"CLDPROP_L2_VIIRS_SNPP\", \"VNP03MOD\", \"viirs_temp\"]] # Suomi NPP\n",
    "sats = [[\"VJ102IMG\",\"VJ103IMG\",\"CLDPROP_L2_VIIRS_NOAA20\", \"VJ103MOD\", \"viirs_temp\"]] # NOAA-2\n",
    "\n",
    "if case2 == \"0313\":\n",
    "    ### Satelite DOY\n",
    "    doy = \"73\"\n",
    "    \n",
    "    ### Times of the satellite overpasses in UTC###\n",
    "    #times = [[\"0212\",\"0354\",\"0530\",\"0712\",\"0854\",\"1030\",\"1212\"], # NOAA-20\n",
    "    #         [\"0118\",\"0300\",\"0442\",\"0624\",\"0800\",\"0942\",\"1124\",\"1300\"]] # Suomi NPP\n",
    "    times = [[\"1212\"]] # NOAA-20\n",
    "    \n",
    "elif case2 == \"0328\":\n",
    "    ### Satelite DOY\n",
    "    doy = \"88\"\n",
    "    \n",
    "    ### Times of the satellite overpasses in UTC###\n",
    "    times = [[\"0412\",\"0554\",\"0730\",\"0912\",\"1048\",\"1230\"], # NOAA-20\n",
    "             [\"0324\",\"0500\",\"0642\",\"0818\",\"1000\",\"1142\",\"1318\",\"1324\"]] # Suomi NPP\n",
    "    \n",
    "else:\n",
    "    ### Satelite DOY\n",
    "    doy = \"89\"\n",
    "    \n",
    "    ### Times of the satellite overpasses in UTC###\n",
    "    times = [[\"0030\",\"0036\",\"0212\",\"0348\",\"0354\",\"0530\",\"0536\",\"0712\",\"0848\",\"0854\",\"1030\",\"1036\"], # NOAA-20\n",
    "             [\"0118\",\"0124\",\"0300\",\"0306\",\"0442\",\"0624\",\"0800\",\"0806\",\"0942\",\"1118\",\"1124\"]] # Suomi NPP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = 9.81\n",
    "\n",
    "ll_lona = np.zeros(2)\n",
    "ll_lata = np.zeros(2)\n",
    "ul_lona = np.zeros(2)\n",
    "ul_lata = np.zeros(2)\n",
    "ur_lona = np.zeros(2)\n",
    "ur_lata = np.zeros(2)\n",
    "lr_lona = np.zeros(2)\n",
    "lr_lata = np.zeros(2)\n",
    "\n",
    "for i in np.arange(2):\n",
    "    if i == 0:\n",
    "        wrf_files = '/glade/campaign/uwyo/wyom0122/lam/new_runs/031220/wrfinput_d01'\n",
    "        wrf_files2 = '/glade/campaign/uwyo/wyom0122/lam/new_runs/032820/wrfinput_d01'\n",
    "    elif i == 1:\n",
    "        wrf_files = '/glade/campaign/uwyo/wyom0122/lam/new_runs/031220/wrfinput_d02'\n",
    "\n",
    "    print ('Reading WRF files...')\n",
    "\n",
    "    s1 = Dataset(wrf_files)\n",
    "\n",
    "    if i == 0:\n",
    "        s2 = Dataset(wrf_files2)\n",
    "\n",
    "    lat2 = getvar(s1,\"lat\",meta=False)\n",
    "    lon2 = getvar(s1,\"lon\",meta=False)\n",
    "    hgt = getvar(s1,\"HGT\",meta=False)\n",
    "    hgt[hgt==0.0] = -100.\n",
    "    if i == 0:\n",
    "    \tlat_d01 = lat2\n",
    "    \tlon_d01 = lon2\n",
    "    \thgt_d01 = hgt\n",
    "    \tseaice = getvar(s1,\"SEAICE\",meta=False)\n",
    "    \tseaice_d01 = seaice\n",
    "    \tseaice2 = getvar(s2,\"SEAICE\",meta=False)\n",
    "    \tseaice2_d01 = seaice2\n",
    "    elif i == 1:\n",
    "    \tlat_d02 = lat2\n",
    "    \tlon_d02 = lon2\n",
    "    \thgt_d02 = hgt\n",
    "\n",
    "    if i < 2:\n",
    "\t    ll_lon = lon2[0,0]\n",
    "\t    ll_lat = lat2[0,0]\n",
    "\t    ul_lon = lon2[-1,0]\n",
    "\t    ul_lat = lat2[-1,0]\n",
    "\t    ur_lon = lon2[-1,-1]\n",
    "\t    ur_lat = lat2[-1,-1]\n",
    "\t    lr_lon = lon2[0,-1]\n",
    "\t    lr_lat = lat2[0,-1]\n",
    "\n",
    "\t    ll_lona[i] = ll_lon\n",
    "\t    ll_lata[i] = ll_lat\n",
    "\t    ul_lona[i] = ul_lon\n",
    "\t    ul_lata[i] = ul_lat\n",
    "\t    ur_lona[i] = ur_lon\n",
    "\t    ur_lata[i] = ur_lat\n",
    "\t    lr_lona[i] = lr_lon\n",
    "\t    lr_lata[i] = lr_lat\n",
    "\n",
    "ll_lon = ll_lona[0]\n",
    "ll_lat = ll_lata[0]\n",
    "ur_lon = ur_lona[0]\n",
    "ur_lat = ur_lata[0]\n",
    "\n",
    "# center point of overall domain\n",
    "ref_lat = (ll_lat + ur_lat)/2.\n",
    "#ref_lon = -(abs(ll_lon) + abs(ur_lon))/2.\n",
    "ref_lon = 8.5\n",
    "\n",
    "# create map subplots\n",
    "# 1000 MB\n",
    "fig = plt.figure(figsize=(20,15))\n",
    "ax = fig.add_axes([0.05,0.05,0.9,0.9])\n",
    "#m = Basemap(llcrnrlon=ll_lon,llcrnrlat=ll_lat,urcrnrlon=ur_lon,urcrnrlat=ur_lat,lon_0=ref_lon,lat_0=ref_lat,projection='lcc',resolution='h')\n",
    "m = Basemap(llcrnrlon=ll_lon,llcrnrlat=ll_lat,urcrnrlon=ur_lon,urcrnrlat=ur_lat,lon_0=ref_lon,lat_0=ref_lat,projection='cass',resolution='h')\n",
    "\n",
    "x_d01, y_d01 = m(lon_d01,lat_d01)\n",
    "x_d02, y_d02 = m(lon_d02,lat_d02)\n",
    "\n",
    "#my_cmap = cmo.deep\n",
    "#my_cmap.set_under('lightskyblue')\n",
    "#my_cmap2 = cmo.thermal\n",
    "\n",
    "#levels = np.arange(0,2300,100)\n",
    "#CS2 = m.contourf(x_d01,y_d01,hgt_d01,levels,\n",
    "#                   cmap=my_cmap,\n",
    "#\t\t   extend='both',\n",
    "#                  alpha=1.0,\n",
    "#                   origin='lower',\n",
    "#                   zorder=1)\n",
    "\n",
    "#CS2 = m.contourf(x_d02,y_d02,hgt_d02,levels,\n",
    "#                   cmap=my_cmap,\n",
    "#                   extend='both',\n",
    "#                   alpha=1.0,\n",
    "#                   origin='lower',\n",
    "#                   zorder=2)\n",
    "\n",
    "###  Loop over all satellites ###\n",
    "sat_i = -1 # index to get satellite\n",
    "for sat in sats:\n",
    "    print ('Doing ', sat)\n",
    "    sat_i += 1 # increase to get next satellite\n",
    "    coord_i = -1 # Index to get coordinate\n",
    "\n",
    "    ### Loop over all times ###\n",
    "    for time in times[sat_i]:\n",
    "        print ('Doing ', time)\n",
    "        coord_i += 1 # increase coordinate index\n",
    "        lon_0, lat_0 = ref_lon, ref_lat # lon and lat at this time\n",
    "\n",
    "        ### all files ###\n",
    "        print (head+sat[0]+\"/\"+sat[0]+\".A20200\" + doy + \".\"+time)\n",
    "        file = sorted(glob.glob(head+sat[0]+\"/\"+sat[0]+\".A20200\" + doy + \".\"+time+\"*\"))\n",
    "        geo_file = sorted(glob.glob(head+sat[1]+\"/\"+sat[1]+\".A20200\" + doy + \".\"+time+\"*\"))\n",
    "        cloud_file = sorted(glob.glob(head+sat[2]+\"/\"+sat[2]+\".A20200\" + doy + \".\"+time+\"*\"))\n",
    "        cloud_geo_file = sorted(glob.glob(head+sat[3]+\"/\"+sat[3]+\".A20200\" + doy + \".\"+time+\"*\"))\n",
    "\n",
    "        ### Import sat files using external script ###\n",
    "        if sat[0][0] == \"V\": # VIIRS\n",
    "            data, lon, lat = import_VNP02.import_VNP02(file[0],geo_file[0]) # Imagery\n",
    "            cth, ctt, cwp, cot, cer, cph, cloud_lon, cloud_lat = import_cldprop_l2.import_cldprop_l2(cloud_file[0], cloud_geo_file[0]) # cloud retrivals\n",
    "            cth_uc, ctt_uc, cwp_uc, cot_uc, cer_uc, cloud_uc_lon, cloud_uc_lat = import_cldprop_uc_l2.import_cldprop_uc_l2(cloud_file[0], cloud_geo_file[0]) # cloud retrivals\n",
    "            \n",
    "        ### Calculate coordinates of corner of LES domain ###\n",
    "        #g = pyproj.Geod(ellps='WGS84')\n",
    "        #diagonal = (2*(frame_size**2))**(0.5)/2 # distance from midpoint of domain to corners\n",
    "        #urc = np.asarray(g.fwd(lon_0, lat_0, 45, diagonal)[0:2]) # upper right corner [lon, lat]\n",
    "        #lrc = np.asarray(g.fwd(lon_0, lat_0, 135, diagonal)[0:2]) # lower right corner [lon, lat]\n",
    "        #llc = np.asarray(g.fwd(lon_0, lat_0, 225, diagonal)[0:2]) # lower left corner [lon, lat]\n",
    "        #ulc = np.asarray(g.fwd(lon_0, lat_0, 315, diagonal)[0:2]) # upper left corner [lon, lat]\n",
    "\n",
    "        ### Mask cloud retrival pixels outside LES domain ###\n",
    "        ### This cloud mask is not perfect especially if the domain is larger ###\n",
    "        ### Because the longitude changes along the sides ###\n",
    "        #mask_cloud = np.logical_or.reduce((cloud_lon>(urc[0]+lrc[0])/2, cloud_lat>urc[1], cloud_lon<(ulc[0]+llc[0])/2, cloud_lat<llc[1]))\n",
    "        mask_cloud = np.logical_or.reduce((cloud_lon>(ur_lon+lr_lon)/2, cloud_lat>ur_lat, cloud_lon<(ul_lon+ll_lon)/2, cloud_lat<ll_lat))\n",
    "        cth[mask_cloud],ctt[mask_cloud],cwp[mask_cloud],cer[mask_cloud],cot[mask_cloud],cph[mask_cloud] = np.nan, np.nan, np.nan, np.nan, np.nan, np.nan\n",
    "        cth_uc[mask_cloud],ctt_uc[mask_cloud],cwp_uc[mask_cloud],cer_uc[mask_cloud],cot_uc[mask_cloud] = np.nan, np.nan, np.nan, np.nan, np.nan\n",
    "\n",
    "### Color bar for infrared imagery\n",
    "vmin = 2\n",
    "vmax = 6\n",
    "color = plt.cm.get_cmap('binary')\n",
    "colorlist  = [color(i) for i in range(color.N)]\n",
    "color = mpl.colors.LinearSegmentedColormap.from_list('Custom cmap', colorlist, color.N)\n",
    "bounds = np.arange(2,6.1,0.1)\n",
    "norm = mpl.colors.BoundaryNorm(bounds, color.N)\n",
    "colormesh = m.pcolormesh(lon, lat, data, cmap=color, latlon=True,vmin=vmin, vmax=vmax,zorder=2)\n",
    "\n",
    "seaice_d01[0:300,:] = np.nan\n",
    "seaice2_d01[0:300,:] = np.nan\n",
    "\n",
    "greenland_ice_del = np.where((lon_d01.ravel()<-19.) & (lat_d01.ravel()>72.))[0]\n",
    "greenland_ice_del2 = np.where((lon_d01.ravel()<-13.) & (lat_d01.ravel()>74.))[0]\n",
    "greenland_ice_del3 = np.where((lon_d01.ravel()<-10.) & (lat_d01.ravel()>80.))[0]\n",
    "seaice_d01_rav = seaice_d01.ravel()\n",
    "seaice2_d01_rav = seaice2_d01.ravel()\n",
    "seaice_d01_rav[greenland_ice_del] = np.nan\n",
    "seaice2_d01_rav[greenland_ice_del] = np.nan\n",
    "seaice_d01_rav[greenland_ice_del2] = np.nan\n",
    "seaice2_d01_rav[greenland_ice_del2] = np.nan\n",
    "seaice_d01_rav[greenland_ice_del3] = np.nan\n",
    "seaice2_d01_rav[greenland_ice_del3] = np.nan\n",
    "seaice_d01_unrav = np.reshape(seaice_d01_rav,(np.shape(seaice_d01)[0],np.shape(seaice_d01)[1]))\n",
    "seaice2_d01_unrav = np.reshape(seaice2_d01_rav,(np.shape(seaice_d01)[0],np.shape(seaice_d01)[1]))\n",
    "\n",
    "m.contour(x_d01,y_d01,to_np(seaice_d01_unrav),levels=[90.],colors=['blue'],linewidths=3.,zorder=3)\n",
    "m.contour(x_d01,y_d01,to_np(seaice2_d01_unrav),levels=[90.],colors=['red'],linewidths=3.,zorder=3)\n",
    "\n",
    "# customize map\n",
    "#m.drawmapboundary(fill_color='aqua')\n",
    "m.drawcoastlines(linewidth=2.0)\n",
    "m.drawparallels(np.arange(60.,94.,4.0),labels=[True,False,False,False],fontsize=28,linewidth=3.0)\n",
    "m.drawmeridians(np.arange(-40.,70.,10.0),labels=[False,False,False,True],fontsize=28,linewidth=3.0)\n",
    "#m.drawstates(linewidth=2.0)\n",
    "m.drawcountries(linewidth=2.0)\n",
    "#m.drawrivers(linewidth=3.0,color='blue')\n",
    "#m.fillcontinents(color='coral',lake_color='aqua')\n",
    "\n",
    "# get x/y coordinates from lat/lon\n",
    "ll_x, ll_y = m(ll_lona, ll_lata)\n",
    "ul_x, ul_y = m(ul_lona, ul_lata)\n",
    "ur_x, ur_y = m(ur_lona, ur_lata)\n",
    "lr_x, lr_y = m(lr_lona, lr_lata)\n",
    "\n",
    "# left\n",
    "m.plot([ll_x[1],ul_x[1]],[ll_y[1],ul_y[1]],linewidth=4,color='white',zorder=4)\n",
    "# top\n",
    "m.plot([ul_x[1],ur_x[1]],[ul_y[1],ur_y[1]],linewidth=4,color='white',zorder=4)\n",
    "# right\n",
    "m.plot([ur_x[1],lr_x[1]],[ur_y[1],lr_y[1]],linewidth=4,color='white',zorder=4)\n",
    "# bottom\n",
    "m.plot([lr_x[1],ll_x[1]],[lr_y[1],ll_y[1]],linewidth=4,color='white',zorder=4)\n",
    "font0 = FontProperties()\n",
    "font = font0.copy()\n",
    "font.set_weight('semibold')\n",
    "x1, y1 = m(40.0,78.8)\n",
    "plt.text(x1,y1,'d01\\ndx=3km',fontsize=28,color='black',bbox=dict(facecolor='darkgray', edgecolor='white', alpha=1.0, boxstyle='round'),zorder=5)\n",
    "##x1, y1 = m(-108.2,42.275)\n",
    "##plt.text(x1,y1,'dx=1000m',fontsize=28,color='white')\n",
    "#x2, y2 = m(28.0,74.0)\n",
    "x2, y2 = m(-4.0,68.1)\n",
    "plt.text(x2,y2,'d02\\ndx=1km',fontsize=28,color='black',bbox=dict(facecolor='darkgray', edgecolor='white', alpha=1.0, boxstyle='round'),zorder=5)\n",
    "##x1, y1 = m(-106.4,40.475)\n",
    "##plt.text(x1,y1,'dx=111m',fontsize=28,color='white')\n",
    "\n",
    "#x3, y3 = m(-121.4,46.3)\n",
    "#plt.text(x3,y3,'Mt. Adams',fontsize=28,color='white')\n",
    "#x4, y4 = m(-122.1,45.3)\n",
    "#plt.text(x4,y4,'Mt. Hood',fontsize=28,color='white')\n",
    "#x4, y4 = m(-121.7,44.7)\n",
    "#plt.text(x4,y4,'Mt. Jefferson',fontsize=28,color='white')\n",
    "\n",
    "#x4, y4 = m(-121.95,45.75)\n",
    "#plt.text(x4,y4,'Columbia',fontsize=28,color='white')\n",
    "#x4, y4 = m(-121.85,45.6)\n",
    "#plt.text(x4,y4,'River',fontsize=28,color='white')\n",
    "\n",
    "#x5, y5 = m(-120.8,46.1)\n",
    "#plt.text(x5,y5,'Washington',fontsize=48,color='white')\n",
    "#x5, y5 = m(-120.4,44.8)\n",
    "#plt.text(x5,y5,'Oregon',fontsize=48,color='white')\n",
    "\n",
    "### Locations\n",
    "#wasco_lat = 45.589999\n",
    "#wasco_lon = -120.672044\n",
    "#gordon_lat = 45.515620\n",
    "#gordon_lon = -120.780007\n",
    "#gv_lat = 45.3620\n",
    "#gv_lon = -120.7456\n",
    "\n",
    "x_anx,y_anx = m(15.874211,69.232286)\n",
    "m.scatter(x_anx,y_anx,marker='D',s=240,facecolor='magenta',edgecolor='k',zorder=8)\n",
    "\n",
    "#x1, y1 = m(wasco_lon,wasco_lat)\n",
    "#plt.scatter(x1,y1,marker='d',s=250,color='magenta')\n",
    "#x1, y1 = m(gordon_lon,gordon_lat)\n",
    "#plt.scatter(x1,y1,marker='o',s=250,color='magenta')\n",
    "##x1, y1 = m(gv_lon,gv_lat)\n",
    "##plt.scatter(x1,y1,marker='x',s=250,color='magenta')\n",
    "\n",
    "#cbarax = fig.add_axes([0.8,0.05,0.04,0.9])\n",
    "#cbar = plt.colorbar(CS2)\n",
    "#cbar.set_label('Elevation (m)',size=36,rotation=270,labelpad=50)\n",
    "#cbar.ax.tick_params(labelsize=28)\n",
    "#cbar.set_ticks([1500,2000,2500,3000,3500,4000])\n",
    "\n",
    "\n",
    "#plt.tight_layout()\n",
    "#plt.subplots_adjust(right=0.775)\n",
    "plt.savefig('domains_d02.png',dpi=250,bbox_inches='tight')\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
